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The overlap Dirac operator at nonzero quark chemical potential involves the computation of the 
sign function of a non-Hermitian matrix. In this talk we present iterative Krylov subspace approx- 
imations, with deflation of critical eigenvalues, which we developed to compute the operator on 
large lattices. We compare the accuracy and efficiency of two alternative approximations based 
on the Amoldi and on the two-sided Lanczos method. The short recurrences used in the latter 
method make it faster and more effective for realistic lattice simulations. 
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1. The overlap Dirac operator at nonzero chemical potential 

Although chiral symmetry can not be implemented exactly in a space-time discretization of 
QCD [1], a lattice version of chiral symmetry can be obtained by a renormalization group blocking 
transformation. The ensuing lattice chiral symmetry is embodied by the Ginsparg-Wilson relation 
[2], which is solved by the overlap Dirac operator proposed by Narayanan and Neuberger [3, 4]. 

Astrophysical objects like neutron stars exhibit an abundance of quarks over anti-quarks, and 
the study of QCD in such a background necessitates the introduction of a quark chemical potential 
jU in the QCD Lagrangian. In this situation chiral symmetry still holds and in Ref. [5] two of the 
present authors proposed an extension of the overlap Dirac operator to nonzero chemical potential 
which still satisfies the Ginsparg-Wilson relation, 

D ov Gu) = l + 7 5 sgn(7 5 D w (At)), (1.1) 

where 75 = 71 727374 with 71 , . . . , 74 the Dirac gamma matrices in Euclidean space, sgn is the matrix 
sign function, and 

3 

DM = 1 - K^iT^ + Tr) - K(e^T+ + e^T 4 ) (1.2) 
i=i 

is the Wilson Dirac operator at nonzero chemical potential [6] with (Ty) yx = (1 ± Yv)U Xt ± v 8 y ^± v , 
K = 1/(8 + 2ra w ), m w £ (—2,0) and U x ^± v £ SU(3). The exponential factors e ±fl implement the 
quark chemical potential on the lattice. 

For /I 7^ the argument y^D^^) of the sign function becomes non-Hermitian, and one is 
faced with the problem of defining and computing the sign function of a non-Hermitian matrix. If 
the matrix A of dimension N is diagonalizable, then A = U diag{A,}?7 _1 with eigenvalues A,- and 
eigenvector matrix U, and a function f(A) can be computed using the spectral definition [7] 

f(A) = Udmg{f(X i )}U- 1 . (1.3) 

If A is not diagonalizable one can still apply an extension of the spectral definition using the Jor- 
dan canonical form. As the eigenvalues of a general matrix A can be complex, the sign function 
computed using Eq. (1.3) requires the sign of a complex number, which can be defined by 

z 

sgn(z) = — = = sgn(Rez) , (1.4) 

VZ Z 

where the branch cut of the square root is chosen along the negative real axis. This definition 
ensures that (sgnz) 2 = 1 and gives the usual result when z is real. It is straightforward to check that 
this definition ensures that (sgnA) 2 = 1. Therefore the overlap Dirac operator at fi / satisfies the 
Ginsparg-Wilson relation, and the operator can have exact zero modes with definite chirality. The 
main properties of the operator were discussed in Ref. [8]. 

Since its introduction the operator has been validated in a number of studies: Its definition is 
consistent with that of domain wall fermions at fi 7^ when the extension of the fifth dimension 
is taken to infinity and its lattice spacing taken to zero [8], its microscopic density and first peak 
computed from quenched lattice simulations on a 4 4 lattice agree with the predictions of non- 
Hermitian chiral random matrix theory [5, 9], and the free fermion energy density was shown to 
have the correct continuum limit [10, 11]. 
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2. Arnoldi approximation for a function of a non-Hermitian matrix 

The exact computation of sgn(A) using the spectral definition (1.3) is only feasible for small 
lattice volumes as the memory requirements to store the full matrix and the computation time 
needed to perform a full diagonalization become prohibitively large for realistic lattice volumes. 
Therefore it is necessary to develop iterative Krylov subspace methods to compute f{A)b for vec- 
tors b. For Hermitian matrices the corresponding methods have already reached a high level of 
sophistication and are widely applied, but for non-Hermitian matrices such methods are novel, ex- 
cept for the special cases of the inverse and the exponential function. The Arnoldi approximation 
discussed below was first introduced in Ref. [12], while the two-sided Lanczos approximation of 
Sec. 4 was presented in Ref. [13]. 

Such iterative methods are based on the observation that the unique polynomial Pk(z) of order 
K < N — 1 which interpolates f(z) at all the eigenvalues of A satisfies the equality 

P K (A)b = f(A)b for any vector b, (2. 1) 

as follows from the definition (1.3). Therefore it is an obvious endeavor to construct a good low- 
degree polynomial approximation for y = f(A)b, taking into account the spectrum of A and the 
decomposition of b in the eigenvectors of A. 

In order to construct such a low order polynomial approximation to f(A)b we consider the 
Krylov subspace K k (A,b) = span(b,Ab,A 2 b, . . . ,A k ~ l b). By definition this subspace contains all 
the vectors resulting from the action of an arbitrary polynomial of degree < k — 1 in A on the source 
vector b. One of these vectors, namely the orthogonal projection of f(A)b on the Krylov subspace, 
will minimize | \P k _\(A)b — f{A)b\ \ over all polynomials of degree <k — \. 

The Arnoldi method uses the recursive scheme 

AV k = V k H k + p k v k+1 el (2.2) 

with 

V^AV k = H k (2.3) 

to build an orthonormal basis V k = (vi, . . . ,v k ) in K k (A,b), where v\ =b/\b\, [} k is the normalization 
of Vfc + i, and e k is the k-th basis vector in C k . The matrix H k is a Hessenberg matrix (upper triangular 
+ one subdiagonal) of dimension k, whose eigenvalues are called the Ritz values of A w.r.t. K k (A,b). 

Once the Arnoldi basis has been constructed, the projection of y = f(A)b on K k (A,b) can be 
written as 

Vproj = V k V}f(A)b = V k Vlf(A)V k vlb . (2.4) 

This formal expression requires the knowledge of the exact solution and is therefore of no practical 
use. However, using the Ritz approximation 

v£f{A)V k ^f(H k ), (2.5) 

which is based on Eq. (2.3), allows us to approximate the projection by 

y V roi~y=\b\V k f(H k )ex. (2.6) 
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By construction this approximation is an element of the Krylov subspace Kk(A,b), and it can be 
shown that the implicit polynomial constructed by this approximation interpolates f(x) at the Ritz 
values of A w.r.t. Kk{A,b) [14]. Note that only the first column of f{H%) is needed in Eq. (2.6). 

This approximation reduces the problem to the computation of with dim/4 dimA, 

which makes it very useful for practical use. The inner sign function /(//*) will then be computed 
with one's method of choice. This could be an exact spectral decomposition if k is small, or some 
suitable approximation method, e.g., Roberts' matrix-iterative method for the sign function, 



which converges quadratically to sgn 

3. Sign function and deflation 

When computing the sign function of the matrix A an additional problem occurs when the ma- 
trix has small eigenvalues, as large Krylov subspaces are then required to achieve a good accuracy. 
The reason is that it is not possible to approximate / well by a low-degree polynomial over the 
entire spectrum of A if it varies too rapidly in a small subregion. A solution to this problem is to 
use the exact value of / for a number of critical eigenvalues of A. Over the remaining part of the 
spectrum / should then behave well enough to allow for a low-degree polynomial approximation. 

Implementing this so-called deflation is straightforward in the Hermitian case, where it is 
based on the fact that any number of eigenvectors span a subspace orthogonal to the remaining 
eigenvectors. In the non-Hermitian case the eigenvectors of A are no longer orthogonal, and a more 
involved approach is needed. We have previously proposed two alternative deflation variants for 
this case [12]. Herein we only present the left-right (LR) deflation and refer to Ref. [12] for the 
details of the Schur deflation. 

The method needs the left and right eigenvectors belonging to the m critical eigenvectors, 



where A m is the diagonal eigenvalue matrix for the m critical eigenvalues and R m = (r\ , . . . , r m ) and 
L m = (£i,... ,£ m ) are the matrices of right and left eigenvectors, respectively. These matrices can 
be made bi-orthonormal, i.e., Z4,/? m = I m , such that R m Lj n is an oblique projector on the subspace 
£l m spanned by the right eigenvectors. 
If the source vector is decomposed as 






(3.1) 



b = b ll +b e 



(3.2) 



where bu = R m Lj n b is an oblique projection of b on Q. m and b e = b — bn, then the matrix function 
can be written as 



f{A)b = f{A)R m V m b + f{A)b e 




(3.3) 



exact 
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The contribution of the first term can be computed exactly, while the second term can be approx- 
imated by the Arnoldi method described in Sec. 2 applied to the Krylov subspace K k (A,b Q ). The 
final approximation is then given by 

f(A)b « R m f{K m )L} n b + \b e \V k f(H k ) ei . (3.4) 

As before the function f{H k ) of the internal matrix should be computed with a suitable method. The 
critical eigenvalues with their left and right eigenvectors should be computed once for all source 
vectors b. 



4. Two-sided Lanczos approximation 



The Arnoldi method described in Sec. 2 suffers from the long recurrences used to orthogonal- 
ize the Arnoldi basis. As an alternative we now consider the two-sided Lanczos method which uses 
short recurrences at the cost of giving up orthogonality for bi-orthogonality. 

Consider the two Krylov subspaces K k (A,b) and K k (A* ,b) and construct two bi-orthogonal 
^V k = I k and the matrix G k = W k \ 



bases V k and W k such that W k V k = I k and the matrix G k = W k AV k is tridiagonal with 



G k = W?AV k 



feci 7i 
Pi a 2 '•• 
'■■ '■■ 



\ 







(4.1) 



(4.2) 



■ ■■ ■■ ■■ Yk-i 
\0 ••• a k J 

It can be shown that V k and W k can be built using the short recurrence relations 

( ftvf+i = (A — cti)vi - Ji-Wi-i , 
1 1 Wf+i = - a*) Wi - ftliW/-! , 

where vi = w\ = b/\b\, a,- = wjAv,-, and /3; and Ji are determined from the normalization condition 
1. 

The matrix V k W k is an oblique projector on K k (A,b), such that an oblique projection of f{A)b 
on K k (A,b) is obtained using 

y ~ Jobi = V k w}f(A)V k w}b . (4.3) 

In analogy to the Arnoldi method we introduce the approximation W, f{A)V k ~ f{G k ) such that 

y u^y=\b\V k f(G k )e l . (4.4) 

The approximation y G K k (A,b), and the problem is now reduced to the computation of f{G k ) with 
dim G,t -C dim A. 

If the matrix A has small eigenvalues, deflation will again be necessary when computing the 
sign function. To implement the LR-deflation in this case one constructs two bi-orthogonal bases V k 
and W k in K k (A, b^ ) and K k (A^,b^), where the directions along R m , respectively L m , have been fully 
deflated from b, i.e., b^ = (1 - R m L] n )b and = (1 - L m R]^)b. With LR-deflation the function 
approximation becomes 

/(A)* « R m f(A m )Ll l b+ \bl\V k f{G k ) ei . (4.5) 
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5. Results 

In this section we compare the results obtained with both methods. In an initial deflation 
phase the right and left eigenvectors of )^D w (^l) corresponding to the eigenvalues with smallest 
magnitude are determined using ARPACK. 1 With these exact eigenvectors the final approximation 
to D ov (n)b of Eq. (1.1) is computed as described in the previous sections. 

In Fig. 1 we compare the convergence rate for both methods, by plotting the accuracy versus 
the Krylov subspace size, and observe that the Arnoldi method has a somewhat better accuracy 
than the two-sided Lanczos method. To achieve the same accuracy, the Krylov subspace of the 
two-sided Lanczos method has to be chosen about 20% larger than in the Arnoldi method. 

However, the speed of the short recunences by far makes up for this as can be seen from 
Fig. 2, where we show the accuracy as a function of the required CPU-time. The two-sided Lanczos 
method is clearly more efficient than the Arnoldi method, and the advantage gained by the short 
recurrences increases as the volume grows. The dotted lines show the time needed to build the 
basis in the Krylov subspace, while the full lines represent the total CPU-time used by the iterative 
method (without the deflation time). The construction of the basis is much faster for the two-sided 
Lanczos method (~ Nk) than in the Arnoldi case (~ Nk 2 ). For the Arnoldi method the time needed 
to construct the Arnoldi basis is dominating, while for the two-sided Lanczos this time can almost 
be neglected compared to that needed to compute the sign function of the inner matrix. Methods to 
boost the computation of the sign function of the inner matrix are the subject of a current study. An 
additional advantage of the short recurrences is their possible implementation with small memory 
footprint if a two-pass procedure is used to compute Eq. (4.5). 

To make the method practical for large-scale lattice simulations it will be important to improve 
the deflation phase, even though first tests on larger lattices (16 3 x 32) with realistic parameter 
values seem to indicate that the number of deflated eigenvalues can be taken much smaller than in 
the 4 4 and 6 4 lattices investigated herein. 

'Typical deflation times on an Intel Core 2 Duo 2.33GHz workstation were t = 27.5s for m = 32 on the 4 4 lattice 
and t = 1713s for m = 128 on the 6 4 lattice. 
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Figure 2: Accuracy versus CPU-time (in seconds), on an Intel Core 2 Duo 2.33GHz workstation, for the 
Arnoldi (red) and the two-sided Lanczos method (blue) for a 4 4 (left) and 6 4 (right) lattice. 

The methods discussed above are also currently being tested to compute eigenvalues of the 
overlap operator. First results are encouraging and clearly show the superiority of the two-sided 
Lanczos method. 
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